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ABSTRACT 

Methods  for  simulating  dependent  sequences  of  contin- 
uous positive-valued  random  variables  with  exponential 
uniform,  Gamma,  and  mixed  exponential  marginal  distributions 
are  given.   In  most  cases  the  sequences  are  first-order, 
linear  autoregressive,  Markovian  processes.   A  very  broad 
two-parameter  family  of  this  type,  GNEAR(l),  with  exponential 
marginals  and  both  positive  and  negative  correlation  is 
defined  and  its  transformation  to  a  similar  multiplicative 
process  with  uniform  marginals  is  given.   It  is  shown  that 
for  a  subclass  of  this  two-parameter  family  extension  to  mixed 
exponential  marginals  is  possible,  giving  a  model  of  broad 
applicability  for  analyzing  data  and  modelling  stochastic  sys- 
tems, although  negative  correlation  is  more  difficult  to  obtain 
than  in  the  exponential  case.   Finally,  two  schemes  for  auto- 
regressive sequences  with  Gamma  distributed  marginals  are  out- 
lined.  Efficient  simulation  of  some  of  these  schemes  is  discussed, 


1.    INTRODUCTION 

In  a  recent  series  of  papers  [1,2,3,4,5,6,7,8,9]  some 
simple  models  have  been  derived  for  stationary  dependent 
sequences  of  positive,  continuous  random  variables  with  given 
first-order  marginal  distributions.   In  general  the  dependency 
structure,  as  measured  by  second-order  joint  moments  (serial 
correlations)  mimics  that  of  the  usual  linear  mixed  auto- 
regressive-moving  average  (ARMA)  models  which  have  been  used 
for  so  long  in  time-series  analysis.   In  the  ARMA  models, 
which  are  defined  quite  generally,  there  is  in  usage  an 
implicit  assumption  of  marginal  normality  of  the  random  vari- 
ables.  This  is  clearly  not  the  case  if  the  random  variables 
are  positive,  say  the  times  between  events  in  a  series  of 
events  (Cox  and  Lewis  [10]}  or  the  successive  response  times 
at  a  computer  terminal.   Thus  the  new  models  are  derived  to 
accommodate  situations  in  which  the  dependent  random  variables 
have,  for  instance,  exponential,  Gamma,  uniform  and  mixed 
exponential  marginal  distributions.   The  exponential  case  is 
the  most  highly  developed,  with  the  nomenclature  (Lawrance  and 
Lewis  [4])  EARMA  (p,q)  (exponential  process  with  mixed  moving 
average-autoregressive  structures  of  orders  p  and  q  respec- 
tively) and  NEARMA(p,q)  (new  EARMA(p,q)).   A  generalization 
to  extend  the  range  of  attainable  autocorrelations  to  nega- 
tive values  has  been  defined  with  the  nomenclature  GNEARMA(p,q) 

The  development  of  the  probabilistic  properties  of 
these  processes  is  given  in  the  referenced  papers,  applications 


to  queueing  models  and  computer  systems  modelling  by  Lewis 
and  Shedler  [11]  and  Jacobs  [12,13],  while  development  of 
estimation  and  testing  procedures  has  just  begun. 

The  object  of  the  present  paper  is  to  define  and 
discuss  the  simulation  of  the  processes  on  digital  computers, 
though  for  the  sake  of  brevity  only  the  first-order  Markovian, 
autoregressive  case  is  considered.   The  simplicity  of  structure 
of  these  models — in  general  they  are  linear  additive  mixtures 
of  random  variables — makes  them  ideal  for  this  purpose.   How- 
ever, stationarity  conditions  are  sometimes  difficult  to  derive 
analytically  and  in  some  cases  it  is  not  simple  to  generate 
the  innovation  random  variables  in  the  processes.   A  striking 
example  of  this  is  the  case  of  the  Gamma  first-order  autore- 
gressive process,  given  in  Section  4A,  for  which  an  efficient 
means  of  simulation  was  reported  by  Lawrance  [7]  for  some  para- 
metric values.   This  procedure  carries  over  into  another  Gamma 
process,  the  Gamma-Beta  process,  which  will  be  discussed  in 
Section  4B. 

In  Section  3  it  is  shown  that  a  simple  transformation 
of  the  exponential  sequences  gives  a  direct  multiplicative 
method  for  generating  dependent  processes  with  uniform  marginals. 
These  could  be  the  basis  in  simulations  for  many  other  types 
of  dependent  sequences. 

Finally  the  NMEAR(l)  process  is  detailed  in  Section  4C; 
this  generates  a  first-order  Markovian  process  with  mixed  expo- 
nential marginals.   It  is  useful  for  simulating  situations  in 
which  the  observed  random  variables  are  correlated  and  over- 
dispersed  relative  to  an  exponentially  distributed  random  variable 


2.    EXPONENTIAL  AUTOREGRESSIVE  MARKOVIAN  SEQUENCES 

We  give  here  three  methods  of  generating  first-order 
autoregressive,  Markovian  sequences  with  exponential  marginal 
distributions.   The  first  two  are  defective  in  terms  of  their 
sample  path  properties  (the  first  more  so  than  the  second) 
while  the  third,  NEAR(l)  and  its  generalization  GNEAR(l) ,  is 
satisfactory  in  this  respect  and  is  a  very  rich  model.   The 
defect  of  the  first  two  sequences  is  also  highlighted  by  the 
simulation  procedures  used;  they  can  be  generated  from  one 
sequence  of  exponential  variables. 

The  word  "autoregression"  in  the  context  of  a  stochas- 
tic sequence   {X  }   is  often  used  rather  vaguely.   In  the 
first  place  linear ,  additive  autoregression  is  usually  implied. 
In  the  second  place  first-order  autoregression  can  mean  that 
in  the  defining  equation  for  X    the  previous  value  enters 
explicitly.   Thirdly,  it  can  mean  that  the  conditional  expec- 
tation of   X  ,  given  X   ,  =  x   - ,  is  an  additive  linear 
n'  ^       n-1    n-1 

function  of   x   , ; 
n-1 


E(x|x1=xn)=a+bx,.  (2.1) 

n1  n-1    n-1  n-1 


The  processes  discussed  in  this  paper  are  autoregressive  in 
the  latter  two  senses  and,  except  in  the  case  of  uniform 
marginals,  are  autoregressive  in  a  linear  additive  way.   They 
are  also  Markovian;  the  Markovian  property  (first-order)  means 
that  the  probability  structure  of  X  ,  X   ,,...  ,  given 

X   ,  =  x    is  independent  of   X   n,  X   -,,.... 

n-1    n         ^  n-2 '   n-3 


2A.   The  Exponential  DAR(l)  Process 

A  very  simple  exponential  autoregressive  Markovian 
sequence  is  generated  by  the  equation  (Jacobs  and  Lewis 
[14,15]) 


X   =VX   ,  +  (1  -  V  )E   ,  (2.2) 

n    n  n-1         n   n  \  «■•«■/ 


where  the  V  's.  n  =  1,2,...  are  i.i. d  with 
n 

P{Vn=l}  =  1  -  P{Vn=0}  =  p   and  E  ,  n  =  1,2,...   are,  as 
throughout  the  paper,  independent  exponential  random  variables 
with  parameter   A   and  independent  of  the   V  's;  that  is 


P{E   <  x}  =  1  -  e"Ax,     x  >  0,     A  >  0 
n  —  — 


(2.3) 


=0        ,     x  <  0. 

For  this  process  the  serial  correlations   p,  =  corr(X  ,X  , .  ) 

r  k         n   n+k 

are 

Pk  =  Pk  (2.4) 

and 

E(X  IX   ,  =  x   -)  =  p,x   ,  +  (1  -  p,)/A.  (2.5) 

n1  n-1     n-1     Kl  n-1         pl  ' 

This  process,  which  was  introduced  to  model  discrete  valued 

variables,  is  not  well  suited  to  modelling  continuous  data 

because  runs  of   X  's   with  the  same  value  can  occur  quite 

n  ^ 

frequently  in  the  sample  paths  of  the  process.   This  happens 


when  X  _,   is  picked  successively  in  (2.2)  ,  rather  than  the 
innovation  E  .   Moreover  the  lengths 
values  are  geometrically  distributed. 


innovation  E  .   Moreover  the  lengths  of  the  runs  of  similar 


2B.   The  Exponential  EAR(.l)  Process 

Another  model  is  derived  from  the  usual  linear  model 

X   =  pX   ,  +  e  (2.6) 

n     n-1    n  v    ' 


in  which  the  i.i.d.  innovation  process   {e  }   is  chosen  so 

n 

that  the   X  's  are  marginally  exponential (A) .   Gaver  and 
Lewis  [1]  show  that  for  this  to  be  true,  one  must  have 
0  <  p  <  1   and 


w.p.  1-p, 

w.p.  p  ,  (2.7) 

where   ^Enh  as  previously,  are  i.i.d.  exponential (A) . 

Again   p,  =  p    and   E(X   X   ,  =  x   ,)  =  p,x   ,  +  (l-p,)A, 
3      k  n1  n-1     n-1      1  n-1       1 

as  at  (2.4)  and  (2.5)  for  the  exponential  DAR(l)  model.   The 
difference  is  in  the  sample  paths;  the  EAR(l)  process  simula- 
tions show  runs  of  geometrically  decreasing   X  's,  but  no  runs 
of  constant  value.   The  geometrically  distributed  runs  occur 
when  only   pX   .   is  picked  in  (2.6). 

The  Markov  property  of  the  two  sequences  implies  that 
if  X0   is  chosen  to  be   E-,  an  exponential ( A)  random  variable 
independent  of   E,,  E2,  ...  ,  then   X,,  X2,.-.  forms  a 
stationary  sequence. 
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Naive  inspection  of  the  defining  equations  (2.2),  (2.6) 
and  (2.7)  suggest  that  to  generate  a  stationary  sequence  of 
length   N,  X,,...,  X., ,  (N+l)   i.i.d.  exponential  deviates  and 
N   uniform  variates  (for  the  selection  process)  are  needed. 
However,  the  sequences  can  be  generated  from  only  one  exponen- 
tial sequence;  this  is  possibly  related  to  the  degeneracy  in 
the  processes.   This  method  uses  the  memoryless  property  of 
exponential ( A)  variables,  namely  that  if   E    is  given  to  be 
greater  than  a  constant   y,  then  E   -  y      is  again  exponential ( A) 

Thus  the  algorithm  for  the  EAR(l)  process  is  to  ini- 
tialize by  setting   X_  =  EQ;  subsequently  set   X   =  pX      if 

E   <  x   =  -£n(l-p)/A;  otherwise  set   X   =  pX   ,  +  (E   -  x  ) . 
n  —   p  n      n—  x      n     p 

This  uses  the  fact  that,  from  (2.3),  P{E   £  x  }  =  p. 

Even  greater  efficiency  can  be  obtained,  though  this 
must  be  qualified  by  considerations  as  to  whether   (i)  the 
X  's  are  to  be  generated  one  at  a  time  or  in  an  array;   (ii) 
a  subroutine  is  available  to  generate  exponential  random 
variables  faster  than  can  be  done  by  taking  logarithms  of 
uniform  deviates,  and   (iii)  the  speed  of  division  in  the 
computer  is  short  compared  to  the  time  needed  for  generation 
of  uniform  deviates. 

The  more  efficient  scheme  recycles  uniform  variables, 
i.e.  if   U   is  given  to  be  between  constants   a   and   b, 
where   0  <  a  <  b  <_  1,  then   (U-a)/(b-a)   is  a  uniform  random 
variable.   (Note  that  its  value  is  not  given,  only  that  it  is 
in   (a,b)  )  .   The  expected  number  of  uniform  deviates  required 


to  generate  an   EAR(l)   process  of  length   N  with  this 
algorithm  is   1  +  (l-p)N,  which  is  less  than  the  number   N 
required  to  generate  an  i.i.d.  exponential (X)  sequence.   Also 
the  expected  number  of  logarithms  is  (l-p)N,  while   N 
comparisons  are  always  needed. 

2C.   The  Exponential  NEAR(l)  Process 

A  broader  two-parameter  exponential  sequence  which 
is  a  first-order  autoregressive,  Markovian  process  and  an 
additive  linear  mixture  of  random  variables  is  given  by 
Lawrance  [7]  and  developed  by  Lawrance  and  Lewis  [5].   Called 
NEAR(l),  the  sequence  is  defined  as 


fex 


n-l 


X   =  e   +  { 
n     n 


w.p 


w.p.    1-a 


n  =  1,2, . . .  ,     (2.8) 


where   0  <    a  <_  1   and   0  <    3  £  1   but   a  =  3  +  1-   Also  the 

selection  process  is  done  independently  for  each   n.   It  can 

be  shown  that  for  the   X    to  be  marginally  exponential ( X) 

the  innovation  variable  e        must  be  generated  from  an   E   by 

n  3  n   J 

the  exponential  mixture 


£    = 

n 


n 


w.p.  6  = 


1-3 


1-  (1-  a)3 


n  =  1,2,... 


(2.9) 


a3 


(l-a)3E     w.p.  1-6  =  t rr^- \ 

^         n      ^  l-(l-a) 


providing   a   and   3   are  not  both  equal  to  one.   When   a  =  0 
or   3=0   the   ^xn^   are  i.i«d.  exponential  variables, 
whereas  when   a  =  1   the  EAR(l)  model  given  at  (2.6)  and  (2.7) 
is  obtained.   In  fact  choosing   a   as  a  function  of   6   in  a 
suitable  way,  e.g.   $  —  a,    gives  an  exponential  model  with  a 
full  positive  range  of  serial  correlation  of  order  one,  since 
it  is  easily  shown  that 

Pk  =  (a3)k  .  (2.10) 

Again 

E(Xn|Xn-1  =  xn_x)  =  a6xn_1  +  Q-a3)/A 

=  P1xn_1  +  (l-p-^/A  2.11) 

and  Xn  =  En   gives  a  stationary  sequence.   Thus  the  correla- 
tions and  regressions  are  the  same  as  for  the  first  two  models. 
However  the  NEAR(l)  process  allows  one  to  model  a  broader  class 
of  exponential  sequences,  as  measured  either  by  sample  path 
behavior  or  higher-order  joint  moments;  see  Lawrance  and  Lewis 
[5]  for  details. 

A  particularly  simple  case  occurs  when   3=1;  this 
model,  called  TEAR(l),  is  very  tractable  analytically  and,  as 
will  be  shown  in  Section  4C,  extends  easily  to  the  case  of 
mixed  exponential  distributions  for  the   X  . 

Note  that  in  the  NEAR(l)  process  the  innovation   e 
is  always  present  unless   a  =  1   and  it  is  therefore  not 


possible  to  simulate  the  stationary  process  with  less  than 
N+l  uniform  variates.   A  detailed  algorithm  is  given  in 
Section  2E  for  a  more  general  case.   Since  for  a  stationary 
array  of  N   X  '  s,  exactly  N+l   uniform  deviates  are  required 
because  of  the  ability  to  recycle  the  uniforms  and  transform 
them  into  exponentials,  it  could  be  advantageous  to  generate 
these  uniform  deviates  in  an  array  which  would  be  replaced 

one  at  a  time  by  the   X  's.   Care  must  be  taken  with  the  re- 

**        n 

cycling  of  the  uniform  variates   U   if  y   =    1-a   is  close  to 
one  or  zero.   In  that  case  it  is  probably  better  for  computa- 
tional reasons  to  use   2 (N+l)   uniform  variates.   Note  that 
a  =  0   gives  the  EAR(l)  process.   When   3=1   a  simpler 
algorithm  can  be  used  since   E    is  no  longer  a  mixture  of 
two  exponentials;   this  important  special  case  is  called  the 
TEAR(l)  exponential  process. 

When   3  =  -^ —  another  one-parameter  sub-class  of  the 

2-a  c 

NEAR(l)  process  is  obtained  with  strikingly  regular  sample 
path  properties.   For  this  process   p^xn  >   xn_i^  =  I/2'*  this 
is  in  striking  contrast  to  the  EAR(l)  process  whose  sample 
paths  show  runs-down  and  the  TEAR(l)  process  whose  sample 
paths  show  runs-up.   This  is  illustrated  in  Figure  1;  all 
sequences  have   p  =  0.75   and  are  transformations  of  the 
same   E    sequence.   The  parameter  space  of  the  NEAR(l)  process 
is  illustrated  in  Figure  2. 
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GNEAR(l)  PROCESS 
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Figure  1.   Sample  paths  of  NEAR(l)  processes,  all  with 
a3  =  p,  =  0.75,  for  different  values  of   a  and   6.   Figure   i 
is  the  EAR(l)  process  (a  =  1,  3  =  0.75),   Figure  B.  is  the 
TEAR(l)  process  (a  =  0.75,  3=1),  Figure  C  is  the  PREAR(l) 
process   a  =  ^~-  =0.857,  and  Figure  D  is  the  REAR(l)  process 


a  =  3  =  (0.75) V 


All  the  sample  paths  are  transformations 


of  the  same  i.i.d.  exponential  sequence 


n  =  0,1, 


,100. 
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/3  =  1 ;  TEAR  ( 1  )  — ► 

irxxxxxxxxxxvxx 


/3  =  0;  I.  I.  d 


Figure  2.   Parameter  space  for  the  NEAR(l)  exponential  auto- 
regressive  process.   The  cases   a  =  0   and/or   3=0   give 
i.i.d.  exponential  sequences  while   3  =  l/(2-a)   gives  the 
partially  reversible  PREAR(l)  process  for  which   P(Xn  <  Xn_^) 
For   a  =  1  we  have  the  original  exponential  process  EAR(l) 
which  tends  to  have  runs-down.   The  TEAR(l)  case,   3=1 
gives  very  simple  analytics  but  exhibits  runs  up.   Also  shown 


is  a  locus  of  constant 


in  this  case 


p,  =  a3  =  .  5 


=  1/2 


12 


2D.   The  Generalized  Exponential  Process  GNEAR(l)  with  Negative 
Correlation 

The  exponential  processes  defined  in  Sections  2A,  2B , 
and  2C  do  not  exhibit  negative  correlation  or  alternation  of 
correlations.   Such  behaviour  is  found  in,  say,  a  normal  linear 
first-order  process  for  which   p.  =  cor(X  ,  X   .)  =  p-1   and 
-1  <  p  <  1   so  that,  for  instance,   p,   can  be  negative.   A 
scheme  for  broadening  the  correlation  structure  of  the  EAR(l) 
process  is  given  in  Gaver  and  Lewis  [1].   However  in  the  ex- 
ponential case  a  much  simpler  alternative  method  is  available. 

Assume  for  simplicity  that   A  =  1.   Now   X   ,   is 

n-1 

a  unit  exponential  variable  and   U   =  F(X   , )  =  1  -  exp(-X   ., ) 

'  n      n— 1  c        n-1 

is  a  uniform  (0,1)  variable,  as  is   1  -  U   =  exp(-X   ,) .   Then 

n     c        n-1 

X^_1  =  F_1(l  -  Un)  =  F~1(l  -  F(Xn_1))  =  -  ln(l    -  expf-X^^)  is 
a  unit  exponential  variable;  in  fact  it  is  the  antithetic 
of   X   ,   which  gives  the  maximum  negative  correlation 
attainable  in  a  bivariate  exponential  distribution: 


r  =  corr(X   .,  X*  ,)  =  1  -  tt2/6  =  -  .6449.  (2.12) 

n-1    n-1  ' 

Now  the  process 


X   = 

n 


or 


e   +3X   n=e   -3^n(l-e  n~)     w.p.a 
n     n-1    n  c 

e  w. p.  1-a  (2 . 13) 

n 


X   =  e   +  U   X   ,  (2.14) 

n    n    n   n-1 
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in  which  the  Un's  are  i.i.d.  with  P{Un=l}  =l-P{Un=0}=a 
gives  a  process  with  autocorrelations  which  alternate  in  sigh. 
In  particular  p,  =  r(a$).  To  combine  this  with  the  positive 
correlation  case  in  a  continuous  way  we  introduce  a  new  para- 
meter p  £[0,1]  and  i.i.d.  indicator  variables  I  ,  indepen- 
dent of  U  ,  from  which  P{l  =1}=1-P{l  =o}=P.  Then 
n  n  n 

the  GNEAR(l)  model  is  defined  as 


-X 

x„   =    £„  +  ur^{l  X   ,  +  (1  -  I  )  [-  ln(l   -  e   n"1)]}         (2.15) 
n    n    n   n  n-±         n 


and  gives  a  complete  range  of  first-order  serial  correlations 


p1   =  a6[p  +  d-p)r]  .  (2.16) 


Higher  lag  correlations  are  more  complicated  and  will  be  given 

elsewhere.   Note,  however,  that   1-p  =  l/(l-r)   gives  a  case 

in  which   X   and   X   ,   are  a  bivariate  exponential  pair  which 
n        n-1  £-         jt 

are  dependent  but  have  zero  correlation. 

Figure  3  give  four  sample  paths  for  the  case   p  =  0 
for  values  of   a   and   6   corresponding  to  those  in  Figure  1, 
which  is  the  case   p  =  1. 

2E.   Algorithms  for  the  GNEAR(l)  Process 

We  give  here  two  algorithms  for  the  GNEAR(l)  process, 
one  based  on  generation  of  uniform  deviates  one  at  a  time, 
the  other  based  on  the  availability  of  subroutines  to  generate 
arrays  of  uniform  and  exponential  variables.   The  sequences 
generated  are  unit  exponentials  (A  =  1).   The  special  case 
a  =  1  (EAR(l))  and   3=1  (TEAR(l))  are  handled  separately  as 
they  will  cause  divides  by  zero  in  the  main  algorithm. 
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Figure  3.   Sample  paths  for  the  GNEAR(l)  process  with   p  =  0, 
corresponding  to  those  in  Figure  1  for  the  NEAR(l)  process 


(p  =  1) .   Here   p. 


-  0.6449  aS. 
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ALGORITHM  GNEAR1A 

This  algorithm  generates  a  sample  of  size  N  from  the  GNEARl 
process.   It  is  based  on  the  generation  of  uniform  random  numbers  one  at 
the  time,  with  recycling  of  these  uniform  random  numbers  for  further  use. 
It  is  assumed  that  a  subroutine  (UNIFORM)  exists  that  generates  raw 
uniforms.   Input  values  are  N,  ALPHA,  BETA  and  P   the  last  three  taking 
values  in  the  closed  interval  [0,1].   However  ALPHA  =  BETA  =  1  is  not  allowed, 


INPUT  N,  ALPHA,  BETA,  P 
CALL  UNIFORM  (U) 


X(0)  +   -loge(U) 

/*  conv 

C  <-   (1- ALPHA)  *BETA 

IF  ALPHA  =  0 

THEN 

D  +   1 

ELSE 

D  ^  (1-BETA)/(1-C) 

/*  Generate  uniform  U*/ 

/*  convert  to  exponential*/ 


END  IF 

DO   I  «-  1  to  N 

CALL  UNIFORM  (U) 

IF  U  <  ALPHA     THEN 


ELSE 


END  IF 

IF  U  £  D 

THEN 

ELSE 

END  IF 

X(I)  *■   Z  +  Y 

DO 

CALL  UNIFORM  (V) 

IF  V  <  P  THEN  Y  <-  BETA*  X(I-l) 

ELSE  Y  ^  -  BETA*  log 

°e 

END  IF 

U  +   U/ALPHA 

Y  +   0 

U  <-  (U-ALPHA)/(1-ALPHA) 

Z  *■  -loge(U/D) 

Z  +■  -   C  *  log  (U-D)/(l-D)) 


l-exp(-X(I-l) 


/*  Recycle  U  */ 


/*  Recycle  U  */ 
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ALGORITHM  GNEARlB 

This  algorithm  generates  a  sample  of  size  N   from  the  GNEAR1 
process.   It  is  based  on  the  generation  of  arrays  of  uniform  and  expo- 
nentially distributed  random  numbers.   It  is  assumed  that  subroutines 
UNIFORM  and  EXPON  exist  that  generate  arrays  of  uniform  random  numbers 
and  exponential  random  numbers  respectively.   For  each  GNEAR1  number 
generated  two  raw  uniform  random  numbers  are  used.   The  first  is  recycled 
in  those  cases  where  a  third  uniform  is  needed  to  make  the  selection  with 
probability  P.   Input  values  are  N,  ALPHA,  BETA  and  P   the  last  three 
taking  values  in  the  closed  interval  [0,1],  except  the  case  ALPHA  =  BETA  =  1 


INPUT  N,  ALPHA,  BETA,  P 
C  «-  (1  -  ALPHA)  *BETA 
IF  ALPHA  =  0  THEN 

ELSE 
CALL  UNIFORM  (U,  2*N) 
CALL  EXPON  (E,  N+l) 
X(0)  +■  E(N+1) 
DO   I  «-  1   to  N 

IF  U(I)  <  ALPHA  THEN 


ELSE 


END  IF 
X(I)  +   Z  +  Y 
END  DO 


D  «-  1 

D  +    (1-BETA)/(1-C) 

/^Generate  array  of  2N  uniforms*/ 

/^Generate  array  of  N+l  exponentials*/ 


V  =  U(I) /ALPHA  /*  Recyle  U  */ 

IF  V  <  P  THEN  Y  «-  BETA*  X(I-l) 


ELSE  Y  <-   -   BETA*  log 


l-exp(-X(I-l)) 


END  IF 

Y  f-  0 


END   IF 

IF     U(I+I)    <   D 

THEN 

Z  «-  E(I) 

ELSE 

Z  +  C*E(I) 
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3.    UNIFORM  MARKOVIAN  SEQUENCES,  NUAR(l) 

It  is  convenient  to  have  dependent  sequences  of  random 
variables  with  marginal  distributions  other  than  exponential. 
Before  discussing  other  solutions  to  Equation  (2.8)  we  show 
that  a  simple  transformation  of  the  NEAR(l)  process  gives  a 
two-parameter  family  of  Markovian  random  variables  with 
uniform  marginal  distributions.   It  is  well-known  that  an 
exponential  transformation  of  a  unit  exponential  random  vari- 
able gives  a  uniformly  distributed  random  variable.   Thus  we 
have  from  (2.8)  and  (2.9)  the  multiplicative  model  for  a 
uniform  Markovian  sequence   ^xn^/  n  =  1/2,...  ,  called 
NUAR(l) , 


where 


X   =  < 
n 


e   X^  , 
n   n-1 


n 


e   = 
n 


U 


n 


u 


(1-a)  3 


n 


w.p.   a 
w.p.  (1-a) 

XI   ~~   -L  /  ^  f    *    *  •    / 


1   ft 

w*p-  Y  =  l-(l-a)3 


w.p.  1-Y  = 


a3 


l-(l-a) 3 


(3.1) 


(3.2) 


(3.3) 


n  =  1,2,... 

for   U  ,  n  =  1,2,...,  i.i.d.  uniformly  distributed,  providing 
that   a   and   3   are  not  both  equal  to  one.   Again  if   Xfl 
is  uniformly  distributed  and  independent  of   U, ,  U2 ,  ...  the 
sequence  is  stationary. 
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The  sequence  is  clearly  quite  simply  extended  to  give 

negative  correlation,  as  in  the  GNEAR(l)  process;  in  fact 

X  -,       in  (3.1)  is  just  replaced  by  (1-X   ,).   Algorithms  are 
n— 1  n— l 

easily  obtained  by  adaptation  of  those  given  in  Section  2E 

for  the  exponential  case.   It  remains  to  find  the  correlation 

structure  and  the  regression  of   X   on   X   , . 

3  n       n-1 

To  do  the  former,  let   X*   be  a  NEAR(l)  sequence  with 

A  =  1,  so  that  the  sequence   X    at  (3.1)  is  given  by 

X   =  exp{-X*}.   Now  the  joint  Laplace-Stielt jes  transform  of  X 
n         n  n 

* 

and  Xn_k,  <J>  *,   #   (s,t)  =  E{exp[-  sX*  -  tX*  .]},  is  given 

n   n-k 
by  Lawrance  and  Lewis  [5].   Setting   s  =  t  =  1   in 

^X*   X*   (s,t)   9ives/ 
n'   n-k 


(J)  *       (1,1)  =  E{exp(-X*)  exp(-X*  L)} 
n'   n-k 


-  E(XnXn_k)  .  (3.4) 


Then  using  the  fact  that  for  a  uniform  random  variable 

E(X)  =  1/2   and   var(X)  =  1/12   we  have,  after  simplification, 


3 


pk  =  corr(Xn,Xn_k)  =  — — -£  J 


ae 


2  +  3   i=l  11  +  (1-a) 3  J 

k  =  1,2, . . .  .  (3.5) 

Note  that  this  is  not  simply  a  geometrically  decaying  corre- 
lation sequence,  as  for  the  NEAR(l)  process.   However,  for 
the  important  special  case  when   3=1  we  get 
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pk 


a 


2-a 


k 

k  =  1,2,...  ,  (3.6) 


and  thus  the  serial  correlations   p.   are  the  kth  power  of 
p,,  which  takes  on  any  value  between  0  and  1.   Thus  we  have 
a  particularly  simple  uniform  Markovian  sequence,  although 
the  sample  paths  will  tend  to  have  runs- up. 

A  similar  analysis  given  in  Lawrance  and  Lewis  [5] 
shows  that 

E'\iVi  - »)  -  7  rn  /a-aiei  {1  -  a  +  m&}  (3-7) 

so  that  the  regression  is  not  linear  for  this  Markov  process 
with  uniform  marginals,  unless   $  =  1. 

This  uniform  sequence  could  form  the  basis,  via  a 
probability  integral  transform,  of  many  other  sequences  with 
given  marginals.   The  parametrization   6  =  ^—  is  a  good 
choice  for  a  one  parameter  model  since  this  case  gives  a  sample 
path  which  is  partially  time-reversible  (see  Lawrance  and 
Lewis  [5]),  with  a  balance  of  runs-up  and  runs-down.   However, 
marginal  transformations  do  not  preserve  correlation  structure, 
as  shown  at  (3.5) ,  and  it  is  therefore  useful  to  see  whether 
sequences  with  marginals  other  than  exponential  can  be  generated 
from  (2.8);  this  requires  finding,  if  possible,  a  suitable 
choice  of  innovation  sequence   e  .   The  result  will  be  a  simple 
process  with  autoregressive  Markovian  structure  and  the  desired 
marginal  distribution. 
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4.   MARKOVIAN  SEQUENCES  WITH  SOME  OTHER  MARGINALS 

Although  an  exponential  distribution  is  a  common  assump- 
tion for  positive  random  variables  met  with  in  problems  in 
operations  research,  it  is  too  narrow  an  assumption  to  encom- 
pass many  real  situations.   Therefore  parametric  distribution 
models  are  invoked  which  include  the  exponential  as  a  special 
case  and  which  allow  for  the  modelling  of  data  which  has  greater 
or  lesser  dispersion  than  exponentially  distributed  data.   Two 
commonly  used  models  are 
(i)   the  Gamma (k, A)  distribution  whose  probability  density 
function  is 


-,/-,   \K— I  — AX 

f(x)  =       r(k^ ,         k  >  0;  A  >  0;  x  >  0,  (4.1) 


where   r(k)   is  the  complete  gamma  function,  and 
(ii)   the  (convex)  mixture  of  exponential  random  variables 


-A-.X  ~^?x 

f(x)     =    TT1A1e  +    (l-TT1)e  ,    0    <    Ax    <    A2; 


x  f>    0,    0    <_  tt1    <_   1    .  (4.2) 


The  Gamma  distribution  has  dispersion,  measured  by  the  coef- 
ficient of  variation   C(X)  =  c(X/E(X)),  which  is  greater  than 
the  exponential  value  of  1  if  k  <  0   and  less  than  1  if  k  >  1. 
The  mixed  exponential  always  has   C(X)  >_  1,  the  equality 
occurring  when  the  special  case  of  an  exponential  random 
variable  with  parameters   A,   or   A~   holds. 
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4 A.   The  Gamma  GAR(l)  process 

The  solution  of  the  standard  first-order  autoregres- 
sive  equation  (2.6)  with  stationary  gamma  marginals  defines 
the  GAR(l)  process.   Using  Laplace-Stielt jes  transforms  with 
(2.6)  shows  that  for   X    to  be  Gamma  (k,A) ,  we  must  have 

cj>£(s)  =  E(e~S£)  =  {p  +  (1-p)  ^4   .  (4.3) 


For   k   integer  this  has  an  explicit  inverse.   For  example, 

2 

for   k  =  2   the  innovation   e   is  zero  with  probability  p  , 

is  exponential ( A)  with  probability   2p(l-p)   and  is  Gamma(2,A) 

2 

with  probability  (1-p)  .   It  is  easy  to  show  in  general  that   e 

is  zero  with  probability   p  ,  so  that  the  "zero  defect"  is  not 
serious  for  large   k.   A  method  of  simulating  a  random  variable 
whose  Laplace-Stielt jes  transform  is  equation  (4.3)  was  derived 
by  Lawrance  [7],  using  the  fact  that  this  sequence  arises  in 
a  particular  type  of  shot  noise  process.   From  this  we  have  the 


Gamma  Innovation  Theorem 

Let  N   be  a  Poisson  random  variable  with  parameter 

6  =  -k  £n(p).   Let   U,  ,  U~,  .  .  .  .  U„   be  uniformly  distributed 

12         N  J 

over   (.0,1)   and  independent.   Let  Y,  ,  .  .  .  ,  Y   be  exponential  (A) 
and  independent.   Then  e      can  be  simulated  using 


£ 


N  U 

VYpmifN>0, 
m=l 


=   0  if     N   =    0    .  (4.4) 
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A  proof  is  not  given  here.   Note  that   e   is  zero  with  Preb- 
le 
ability   exp{-k  £n(p)}  =  p  .   Also  the  Poisson  number   N 

of  uniform  and  exponential  random  variables  which  must  be 

generated  for  each  e      has  expected  value   0  =  -k  £n(p) . 

This  will  be  prohibitively  large,  and  the  simulation  will  be 

very  inefficient,  if   k   is  large  and/or   p   is  close  to  zero. 

Neither  of  these  cases  is  serious,  however.   If   k   is  large, 

say  greater  than  50,  the  sequence  is  almost  normal  and  the 

usual  normally  distributed,  AR(1)  linear  process  can  be  used. 

If   p   is  as  small  as  0.001  then   E (N)   is  only   k  x  (6.9078) 

which  is  still  reasonable.   However,  for   p   this  small  the 

sequence  is  approximately  i.i.d.  Gamma  and  acceptance-rejection 

techniques  for  simulating  Gamma  variables  are  known. 

It  is  quite  simple  to  write  algorithms  for  the  GAR(l) 
case  analogous  to  those  in  Section  2E.   It  would  pay  to  have  a 
built-in  routine  for  generating  the  Poisson  variable  which  will 
bypass  further  calculations  if   N  =  0.   In  other  words  routines 
for  generating  Poisson  variates  which  start  by  searching  at  the 
median  of  a  table  of  cumulative  Poisson  probabilities  will  be 
inefficient. 

Unfortunately  the  NEAR(l)  process  does  not  appear  to 

extend  to  the  Gamma  case;  it  can  be  shown  explicitly  that 

there  is  no  innovation   e    in  equation  (2.8)  which  will  make 

n 

X   have  a  Gamma  distribution  with  k  =  2  if  a  =f  1. 

n  ' 

There  is,  however,  another  model,  the  Gamma-Beta  model, 
inspired  by  an  example  in  Verwaart  [16]  and  discovered  inde- 
pendently by  Fishman  [17],  which  is  quite  broad  and  which  can 
be  simulated  using  the  Gamma  Innovation  Theorem. 
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4B.   The  Gamma-Beta  Model,  GBAR(l) 

This  model  is  a  linear  autoregressive  process  with 
random  coefficients  which  includes  the  GAR(l)  process  but  is 
of  limited  practical  use  in  data  analysis  because  its  likeli- 
hood function  is  analytically  untractable.   Nevertheless  it 
could  be  useful  in  simulations,  particularly  if  the  zero-defect 
in  the  GAR(l)  process  is  unacceptable. 

Thus,  we  define  the  stationary  sequence  as 

Xn   =    3BnXn_1   +en  0<e<l,       n=0,    +1,    +2,...,(4.5) 

where  the   B's   are  i.i.d.  and  independent  of   X   ,  ,  and 

n  e  n-1' 

B   ,  for   X    to  have  a  Gamma  (k,A)  distribution,  has  a 
Beta(k-b,b)  distribution  with  mean  E (B  )  =  (k-b)/k.   The  density  is 

/   \  r (k)         k-b-1   ,  ,    \b-l  n      „  ^   t  i a      r\ 

fB  (x)  =  r(k-b)r(b)  x     (1"x)   •     Q  <  x  <  l  ,     (4.6) 

n 

0  <  b  _<  k  . 

The  distribution  of  the  i.i.d.   e's   to  make   X   Gamma  (k,A) 

n  n  ' 

is  still  to  be  determined.   To  do  this  and  to  see  the  rationale 

behind  the  model,  it  is  simplest  to  obtain  the  distribution  of 

B  X   ,  . 
n  n-1 

Now  B    can  be  generated  as   Z,/(Z,+Z2)   where 

Z,   and  Z2   are  independent  and  Gamma  (k-b,A)  and  Gamma  lb, A) 

respectively.   Moreover,   B    is  independent  of  (Z,+Z~)  ,  which 

could  be  used  to  generate  X   n .   Then, 

^  n-1 

BnXn-l   =    {zi/(zi+z2^    ^zi+z2^    =   Zl      is    Gamma    (k-b,A)     •      Using 
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this  fact,  which  can  also  be  shown  analytically,  and  the 
defining  equation  (4.3)  we  have 

E(e     /  =  ^X  (s)  =  ^B  X(3s)<})e  (s)  '  (4'7) 

^  n         n  n-1    n 

so  that  on  using  the  fact  that  for  a  Gamma (k, A)  variable  the 
Laplace-Steiltjes  transform  is  (A/A+s)  ,  we  have 


4>   (s)  =  <j>   (s)/<|>B      (6s) 
n        n      n  n-1 


/    A    \k/A+3slk"b  ,A    „, 

=   (ATi)    (aTs-J  (4.8) 

/    A    \br  i   lk~h                                   (4-9> 

=   (xTi)    [3  +    M-3)  xTi]        • 


Thus   e    is  generated  as  the  sum  of  a  Gamma  (b,A)  variate 
and,  from  (4.3) ,  a  Gamma  innovation  variable  which  can  be 
generated  from  (4.4). 
For  this  model 

p.  =  corr(XnXn+j)  =  (s  t^)3      j  =  0,1,2,... 

and 

E(X^  X„  -  =  x„  n)  =  P,xn  ,  +  (l-p.)E(X  ). 
n   n-1     n-1      1  n-1       1    n 

Also  when   b  =  0   we  have  the  GAR(l)  model.   If   3=1   the 
model  requires  only  a  Gamma  (b,A)  variate  for   e    and  a 


(4.10) 
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Beta(k-b,b)  variate  for  B    in  the  simulation.   The  case 
k  =  1   gives  an  exponential  process  which  is  not  the  same  as 
the  NEAR(l)  process. 

It  should  be  noted  that  while  the  Gamma  Innovation 
Theorem  makes  the  Gamma-Beta  model  tractable  from  a  simulation 
viewpoint,  it  is  still  difficult  to  unite  down  a  likelihood 
function  or  to  get  negative  correlation.   Thus  further 
developments  are  needed  from  the  Gamma  case.  An  algorithm  is 
given  for  the  Gamma  Beta  Model  on  the  next  page. 

4C.   Mixed  Exponential  Markovian  Processes  MEAR(l)  and  TMEAR(l) 

In  addition  to  Gamma  processes,  first-order  autore- 
gressive  Markovian  processes  with  mixed  exponential  marginal 
distributions  can  be  obtained  from  equations  (2.8)  and  (2.9) 
in  two  special  cases,  and  these  sequences  should  be  widely 
useful  in  modelling  stochastic  systems. 
(i)   The  case   a  =  1;  MEAR(l) . 

In  Gaver  and  Lewis  [1J  it  is  shown  that  the  solution 

to  the  Laplace  transform  of   e    for  the  linear  model  (2.6) 

c  n 

is  a  constant   p   plus  a  (generally)  non-convex  mixture  of 
three  exponential  functions.   This  can  be  shown  to  be  a  proper 
density  function  if   p  <_   A,/A2/  but  it  can  also  be  shown  that 
it  is  not  a  density  function  for  all   p   less  than  one  and 
greater  than  or  equal  to  zero.   More  particularly,  Lawrance  [6] 
showed  that  unless   A,   is  much  smaller  than   A2   (and  thus 
the   X   are  very  over-dispersed  relative  to  an  exponential 
random  variable)  a  solution  exists  for   e    for  all   p.   Thus 
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ALGORITHM  GBAR1  (GAMMA-BETA  MODEL) 

This  algorithm  generates  a  sample  of  size  N   from  the  GBAR1 
process  using  array  generators  POISSON,  UNIFORM,  EXPONENTIAL  and  GAMMA 
to  produce  Poisson,  Real  Uniform  (0,1),  Exponential  and  Gamma  distributed 
random  numbers. 

The  following  restrictions  exist  on  the  input  parameters: 

0  1   B  1   K 
0  <  K 

0  <  BETA  <  1 


INPUT  N,  BETA,  B,  K 
TH  +■  -(K-B)  *  log  (BETA) 
CALL  POISSON (PSN,N,TH) 
CALL   GAMMA(Z1,N,K-B) 
CALL   GAMMA (Z2,N,B) 
CALL   GAMMA (X(0),  1,K) 
DO   I  «-  1   to  N 

CALL  UNIFORM (U,PSN (I) 
CALL  EXPON(E,PSN(I)) 
Y  +   0 


/*  Initialization  */ 
/*Generate  N  poisson  deviates  with  parameter  TH  */ 
/^Generate  N  Gammas  with  parameter  K-B  */ 
/^Generate  N  Gammas  with  parameter  B  */ 
/*  Initialize  X(0)   */ 

/*  Generate  PSN(I)  real  (0,1)  uniforms  */ 
/*  Generate  PSN(I)  Unit  Exponentials  */ 


DO  J  +   1  TO  PSN(I) 

Y  *■  Y  +  E(J)  *  BETA  **  U(J) 
END  DO 

BN  «-  Z1(I)/(Z1(I)  +  Z2(I))       /*  Compute  BETA  Deviate  */ 
X(I)  «-  BETA  *  BN  *  X(I-l)  +  Y 
END  DO 
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we  have  a  useful  process,  although  again  the  zero-defect  of 
order   p   is  a  problem.   However,  one  case  which  cannot  be 
simulated  this  way  occurs  when   1/A„  =  0  ,  i.e.   X   is  zero 
with  probability   1-tt,  ,  and  exponential  (A,)  otherwise.   This 
kind  of  situation  occurs  in  practice  as,  e.g.,  the  waiting  time 
for  an  item  in  an  inventory  system.   Fortunately  it  can  be 
handled  in  the  next  case. 

(ii)   The  case   6  =  1;  TMEAR(l). 

When   3=1   in  equation  (2.8) ,  a  mixed  exponential 
process  TMEAR(l)  is  obtained  which  is  extremely  simple  to 

simulate  since  the  innovation   e    is  just  the  mixture  of  two 

n      J 

exponentials  for  all   0  <_   p  <  1.   Moreover,  the  process  has 
no  zero-defect.   As  discussed  above,  the  sample  paths  will 
tend  to  "run  up,"  but  this  is  no  great  problem  unless   p   is 
fairly  large.   Thus  we  have  the  following  Theorem  (Lawrance  and 
Lewis  [18])  which  we  state  without  proof: 

TMEAR(l)  Theorem 

Let  the  first-order  autoregressive,  Markovian  sequence 

{X  }   be  defined  by 

n  J 


X   =  e   +VX   ,,        n=  1,2,3,... 
n    n    n  n-1  '  '  '  ' 


where  for  the  i.i.d.  sequence   V  ,  n  =  1,2,3,..., 

P{V  =1}  =  1  -  P{V  =0 }  =  a   for   0  <  a  <  1.   Then  the  sequence 
n  n  —  n 

{X  }   is  stationary  and  has  a  (convex)  mixed  exponential 
marginal  distribution  with  probability  density  function 
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-X1x  ~^2X 

fx(x)    =    7T1A1e  +    tt2  ,  x    >    0  (4b11) 


where         0    <    \±    <    A,,;    0    <    ^    <    1 ;    u±    +   ^    =    i    ^    =     1_  and    ^    =   _1_      ^ 

if   e    is  i.i.d.  and  has  a  (convex)  mixed  exponential 
distribution  given  by 


-YjX  -y9x 

f£(x)    =    n^e  +   n2Y2e  x  -  ° 

with  y1    >    y2    >    °'       n-^f    n2    >    0,     n1    +    n2    =    1  (4.12) 

where 

\i    =   E(X)     =    tt2    y2    +    TT1y1    ;  (4.13) 

b   =  ux    +    y2    -    ay    ;  (4.14) 

3   =  ux    +   y2    -    y    ;  (4.15) 

a   =  (1-a)  y.^    ;  (4. 16) 

Y1/Y2    =  j{b+/[b2-4a]}    ;  (4.17) 

YQ    =  Tr2yx   +    7rxy2    •  (4>18) 

nl    =    (Y1"Y0)/(Y1~Y2)     '*     n2    =    (Y2~Y0)/(Y2"Y1)  (4.19) 

and   X_   is  independent  of   e, ,  £-,     •••   and  has  probability 
density  function  (4.11). 

Note  that  the  special  cases  where  tt.  =  0   or   tt2  =  1 
give  NEAR(l)  exponential  processes  with  parameters   A~   and 


29 


A,   respectively.   Thus  they  should  be  handled  by  Algorithm  2 
since  they  will  cause  computational  problems.   The  case 
X      =  A9   also  gives  a  NEAR(l)  process  and  is  excluded  for 
similar  reasons.   This  guarantees  that   Yi  >  Y2  "   The 
algorithm  on  the  following  page  implements  this  theorem. 

5.  GENERALIZATIONS 

In  all  of  the  processes  discussed  here  except  the 
exponential  GNEAR(l)  the  correlations  are  non-negative  and 
geometrically  decreasing.   A  particular  scheme  for  obtaining 
negative  correlation  is  given  for  the  exponential  case. 
Another  scheme  for  obtaining  alternating  correlations  which 
are  possibly  negative  and  which  is  broadly  applicable  is 
given  in  Gaver  and  Lewis  [1]  and  in  Lawrance  and  Lewis  [5]. 
Another  problem  is  that  different  types  of  dependence  and 
higher-order  Markovian  dependence  might  be  encountered  in 
data.   Schemes  for  obtaining  mixed  autoregressive  moving 
average  exponential  sequences  where  the  autoregression  has 
order   p   and  the  moving  average  has  order   q   are  given  in 
Lawrance  and  Lewis  [4],   The  mixed  exponential  process 
TMEAR(l)  is  easily  extended  to  give  a  process  with  this  type 
of  extended  correlation  structure.   This  will  be  discussed 
elsewhere. 
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ALGORITHM  TMEAR1 


This  algorithm  generates  a  sample  of  size  N   from  the  TMEAR1 
process.   It  is  based  on  the  generation  of  arrays  of  uniform  and  unit 
exponential  numbers.   Subroutines  UNIFORM  and  EXPON  are  assumed  to  exist 
to  generate  such  arrays.   For  each  TMEAR1  number  generated  two  raw  uniform 
random  numbers  and  one  exponential  number  are  needed. 


INPUT  N,  PI1,  PI2,  MU1,  MU2,  ALPHA 

MU  4-  PI1  *  MU2  +  PI2  *  MU2 
B  «-  MU1  +  MU2  -  ALPHA  *  MU 
A  «*-  (1  -  ALPHA)  *  MU1  *  MU2 
T  +■  SQRT(B  *  B  -  4  *  A) 

G2  <-   .5  *  (B  -  T) 

Gl  +   .5  *  (B  +  T) 

GO  +   PI2  *  MU1  +  PI1  *  MU2 

El  <-    (Gl  -  G0)/G1  -  G2) 

Rl  =  1/G1 

R2  =  1/G2 

CALL  UNIF0RM(U,2*N+1) 

CALL  EXPON (E,N+1) 

IF  U(2N+1)  <  PI1 


THEN 
ELSE 


END  IF 

DO   I  -«-  1   to  N 
IF  U(I)  <  El 


THEN 
ELSE 


END  IF 

IF   U(I+I)  <   ALPHA  THEN 

ELSE 
END  IF 
X(I)  *■  Y  +  Z 
END  DO 


/*  Initialization  */ 


/*  Initialization  */ 
/*  Generate  2N+1  uniforms  */ 
/*  Generate  N+l  unit  exponentials  */ 
X(0)  <-   MU1  *  E(N+1) 
X(0)  *-  MU2  *  E(N+1) 


Z  <-   Rl  *E(I) 
Z  ■«-  R2  *E(I) 

Y  *-  (X(I-l) 

Y  4-  0 
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